%Create a fresh environment
close all
clear all
clc

Tprec = 1/1000;             % Precision used to represent Continous time
t = linspace(-4,4,101);     % Time vector
x = sinc(t);                % Input Magnitude Vector
L = length(x);              % Length of Mag. Vector
c = x.*hann(L)';             % Windowed Mag. Vector
nr = L-1;                   % Number of roots to find

% Seach Precision
N = 8000;
rp = 1/500;
k = 2:N+1;                  % Grid Increment
r = 1.5;                    % Starting rdius

%Initialize Variables
theta = 0; indx = 0; lx = 0; lf_rts=[];zzz=[];

%Compute Initial FFT's
crA= [ 1 fftshift(abs(fft((r+rp).^(0:L-1).*c,N))) 1];
crB= [ 1 fftshift(abs(fft(r.^(0:L-1).*c,N))) 1];

figure();
while (r > 2*rp)
    crC=[1 fftshift(abs(fft((r-rp).^(0:L-1).*c,N))) 1];
    
    % Check for local minima (root)
    Bk = crB(k);
    locs = (Bk < crA(k-1)) & (Bk < crA(k)) & (Bk < crA(k+1)) ...
                & (Bk < crC(k-1)) & (Bk < crC(k)) & (Bk < crC(k+1)) ...
                    & (Bk < crB(k-1)) & (Bk < crB(k+1));
    
    indx = find(locs);              % indices of local minima
    lx = length(indx);              % number of minima found
    
    if (lx)                         % Check if roots found
        theta = -2*pi/N*(indx-1);   % Convert indices into complex numbers
        rts = r*exp(1i*theta);   
        for i=1:length(rts)
            Vroots(i)=laguerre(c,rts(i),1E-14,10000);
        end
        
    end
    
    
    
    %Adjust for next seach
    r = r-rp;
    crA = crB;
    crB = crC;
end
    % Refine roots

%     X = [1 2 3 3 4 5 2]';
Vroots = sort(zzz,1)';
idx = [false;diff(Vroots)<(10^-14)];  %Duplicate precision
Vroots(idx) = [];
%Deflate and solve remaing roots
% if (length(zzz) < nr)
%     pd = roots(deconv(c,poly(zzz)));
%     zzz = [zzz pd'];
% end
% 
% %Plot results
zplane(c);
hold on;
plot(zzz,'*r');
% plot(-1*lf_rts,'*k');
% plot(lf_rts,'*k');
legend('MATLAB roots','LF Roots');
grid on
% 
% figure()
% subplot(3,1,1)
% plot(c);
% subplot(3,1,2)
% plot(0:L-1,-poly(roots(c)))
% subplot(3,1,3)
% plot(0:99,poly(lf_rts))
% 
